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PakalrA 3D Model to Solve the Radiative Transfer Equation 
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ABSTRACT 

We present a new numerical model called "PAKAL" intended to solve the 
radiative transfer equation in a three dimensional (3D) geometry, using the ap- 
proximation for a locally plane parallel atmosphere. Pakal uses pre-calculated 
radial profiles of density and temperature (based on hydrostatic, hydrodynamic 
or MHD models) to compute the emission from 3D source structures with high 
spacial resolution. Then, Pakal solves the radiative transfer equation in a set of 
(3D) ray-paths, going from the source to the observer. Pakal uses a new algo- 
rithm to compute the radiative transfer equation by using an Intelligent System 
consisting of three structures: a cellular automaton; an expert system; and a 
program coordinator. The code outputs can be either two dimensional maps or 
one dimensional profiles, which reproduce the observations with high accuracy, 
giving in this way, detailed physical information about the environment where 
the radiation was generated and/or transmitted. We present the model applied 
to a 3D solar radial geometry, assuming a locally plane-parallel atmosphere, and 
thermal free-free radio emission from a Hydrogen-Helium gas in thermodynamic 
equilibrium. We also present the convergences test of the code. We computed the 
synthetic spectrum of the centimetric - millimetric solar emission and found bet- 
ter agreement with observations (up to 10 4 K at 20 GHz) than previous models 
reported in literature. The stability and convergence test show the high accuracy 
of the code. Finally, Pakal can improve the integration time by up to an order 
of magnitude compared against linear integration codes. 
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1. Introduction 

The observation and study of radio emissions coming from distant sources is a valu- 
able tool to investigate these objects and the medium between them and the observer. For 
instance, by assuming an emission mechanism, we are able to obtain detailed physical prop- 
erties of the observed object as the density, temperature, magnetic field, etc. 

Generally, due to observational limitations, we obtain two dimensional projections in 
the plane of the sky of the emission and/or absorption in the ray paths of each point of 
the observed region inside the telescope field of view. These maps represent not only the 
emitting object, but all the possible flux changes due to the emission and/or absorption 
that may take place in the medium between the source and the observer. Therefore to get 
reliable information through the study of radio emissions, it is necessary to take into account 
the detailed three dimensional (3D) structure of both, the source and the medium. In this 
work we present "Pakal" , a numerical model intended to solve the radiative transfer equa- 
tion, designed to study astronomical objects specially in the millimeter and submillimeter 
wavelengths. 

Now a days, there are few dozens of codes to solve the radiative transfer equation, 
though, each code is designed to solve a very specific problem. As instanc e, there are codes 



to so l ve the radiatiye tran s fer equation in Earth-lik e atmospheres (e. g. lOreopoulos et al. 



2006; Cahalan et al. 2005; Davis fe Cahalan 200 



). More specifica 



ly, the I3RC Monte 

Carlo community m odel of 3D radiative transfer (iCahalan et al.l 120051); the A RTS package 
( Buehler et aHl2005h: Ba ttaglia-Manto vani model faattaglia fe MantovanJ2005h : GRIMALDI 
(IScheirer fe MackdboOlh : MCARaTS Jlwabuchill2006h : SHDOM jEvanslll998f )~and SHARM- 
3D (jLyapustinll2002l ) . are designed to study the dispersion of tele-communication radio waves 
in the Earth atmosphere. These models simulate layers of a plane-parallel atmosphere and 
are based mainly in Monte Carlo techniques. 

In the astrophysics community, there are mainly two branches of codes to simulate the 
emission of stellar atmospheres: 

• Codes to simulate the atmosphere structure (the variation with height of physical 
parameters as density, temperature, etc). 

• Codes to compute the synthetic spectrum. 



Commonly, the codes for stella r atmosphere s simulatio n deal with specific ph ysical 
conditions. For example, ATLAS 1 2 ( ]Kuruczlll979l ); MARCS (IGustafsson et al.lll975l ); and 
PHOENIX (jliauschildt et al" 1999 ) are general propose codes for stellar atmospheres which 
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take into account o nly the emission s from the stellar photospheres. The PANDORA (IVernazza et al 
19761 ) and MULTI (lCarlssonlll992l ) codes simulate stellar atmospheres using conditions sim- 
i lar to the solar atmosphere but only in the region below the corona. Whereas CHANTI 
( IDere et al.l 119971 ) simulates atmospheres with coronal conditions. The chromospheric and 
coronal codes are oriented to reproduce the ultraviolet (U V) and Visible sp e ctrum, and 



therefore fail to reproduc e observations in the radio range (IZirin et al.l Il99ll ; lEwell et al. 
19931 : ISelhorst et al.ll2005h . 



Examples of codes in the second branch (synthe tic radiative spect rum) are: SYNTHE 



flKuruczlll979f ): SPECTRUM flHubenv fc Lanzlll995h : and FANTOM flCavrel et al.lll99lf ). 
These codes are complementary to codes in the former branch and are necessary to compute 
the final stellar spectrum. We note that PANDORA and CHANTI, from the first group, are 
also able to compute the spectrum. 

Some codes are int ended to compute particular stellar atmospheres e. g.: STERNE3 



fBe hara & Jefferv!l2006[ ) for Hydrogen-Deficient Stars; LINE-BY-LI NE METHOD f Shulvak et al 



2004) for stars in early and inte rmediate stage; PRQ2 (jWernerlll986l ) and TLUSTY (IHubeny fc Lanz 
995 ) for hot stars; WM -basic ( iPauldrach et al.ll200ll ) for expandi ng atmospheres; CMFGEN 



flHillier fc Milled Il998f ) for Wolf-Rayet stars; and FASTWIND f lSantolava-Rev et al.lll997f ) 



for stars with high mass loss. 

Pakal is a completely new code which can be applied to any geometry, radiation and 
absorption mechanisms (focused in this work to millimeter and submillimeter wavelengths, 
but easily configurable for other wavelengths). This flexibility is achieved by the means of 
four completely independent modules: the numerical model (Section [2]); geometry model 
(Section l3.2j) ; numerical methods; and physical functions (Section l3.3p . 

Pakal uses a new method to compute the radiative transfer equation in a set of 3D 
ray paths, this is an intelligent system called "Tulum" (Section I3.ip which helps to reduce 
the integration time up to one order of magnitude as compared against direct integration 
codes (Section H]). Pakal is able to compute the contributions to the opacity function of each 
chemical element and its ionization states. To accomplish this, the code needs as input, 
detailed profiles of electron temperature and ion densities (Appendix IA~|) . 



2. Radiative Transfer Theory 



The specific intensity is the most basic entity in radiative transfer theory, and is defined 
as the amount of energy dE passing though an area dA, during a time dt, coming inside a 
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solid angle dw, in an interval of frequency dv, with a direction given by f (jRohlfall986l ) 

dE 



dl u 



dA dt dw dv f ■ n 

where f and h are the direction and normal (to dA) unitary vectors, respectively and can be 
written as, 

f • n = cos 9 — fi, 

where 9 is the angle between f and n. When radiation interacts with matter, crossing a 
distance ds, the change in the specific intensity dl u is equal to the e mission of the medium , 
e v , minus the radiative energy absorbed by the medium, k v I v , this is ( 1Chandrasekharlll960l ): 



dh 
ds 



X 



where k u is the opacity function which depends on the physical properties of the medium. 
Assuming a plane-parallel atmosphere, it is possible to write ds, in terms of the geometric 
distance dx (see Figure H]), 

dx = ds cos(6 ) ) = /ids; 

then, using the optical depth, dr u = —K v dx; and Kirchhoff's law (e„ = k v S u ), Eq. [CJmay be 
written as: 



dl v l v 



St, 



dr v fi fi 

The solution, in the [ri jjy , T2 jU \ optical depth interval (where r^ u > r 2j ^) is 



-{Tl,v—T2,v)/H 



S u {T u )e- {T »- T ^ )/lM d 



(2) 



For solar con ditions, the scattering is negligible in the millimeter and submillimeter 
wavelength range (IVernazza et al.l Il976l ) . Therefore for our purposes Eq. [2] is completely 
valid. 

Assuming < T\\ /j, = 1; and a source function constant in each cell "z" (0 < i < n), 
we integrate Eq. [2] in an array of n consecutive cells (see Figured]), using: 



I„(L<+i) = I u (Li)exp 
+S v {Li + 0.5rfL) ( 1 - exp 



—^(k v {Li) + K(L i+1 )) 



~(K(Li) + K(L i+1 )) 



(3) 



where I (LA is the specific intensity (coming) from the cell "i — 1"; I(L i+ i) is the specific 
intensity (getting out) of cell "i" ; and dL is the integration step. The computation of each 
Li is done by the geometry module (see Section 13721) . 
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3. Pakal Model 

Pakal (the name of the king of Pa lenque in the Mayan Culture) is written in C language 



using an object oriented technique (jSchildtl 119871 ) which allow us to encapsulate sets of 



common properties or functions in libraries. The code is based on four independent modules: 
i) the numerical model, ii) the geometry, iii) numerical methods and iv) physical functions. 

Once the geometry is defined, Pakal generates a series of independent ray paths, from the 
source to the observer, reads pre-defined temperature and density profiles and, if necessary, 
performs an interpolation of the readed values, covering a larger number of points in altitude. 
Then, using an intelligent system called Tulum solves the radiative transfer equation (the 
related algorithms are part of the numerical module). 



3.1. Tulum: The Intelligent System 

The Intelligent System used in Pakal is called Tulum and helps to solve the radiative 
transfer equation in a new and very efficient way. In Figure [2] a schematic diagram of the 
automaton is presented. 

Tulum is formed by three independent components: 

• A coordinator which controls each step of the integration process. The coordinator uses 
the recommendations of the expert system and the states of the cellular automaton to 
decide the next stage of the integration process. 

• An expert system who recommends, based on the current status (position and physical 
conditions), whether or not it is necessary to integrate in this point and, if necessary, 
recommends a change of the integration step size. 

• A cellular automaton, with a set of previously established states, which is able to save 
the current status of the integration process. 

Tulum can integrate numerically any given function (not only the radiative transfer 
equation). The integration process is carried out in the following way: 

1. When the coordinator program receives the spatial coordinates of two contiguous inte- 
gration points (from the 3D geometry module), he looks for the physical conditions at 
these points (from the temperature and density radial models). If necessary, the nu- 
merical module of Pakal automatically interpolates the radial temperature and density 
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models, at the specific points, using either of two classical methods: linear or cubic 
spline interpolations. In this work we use linear interpolation (the cubic spline inter- 
polation fails because the temperature profile has a very large gradient in the Solar 
Transition Zone). 

Once the coordinator knows the physical conditions of the medium, he asks for a rec- 
ommendation to the expert system, and also asks for the present state of the cellular 
automaton. Based on this information, the coordinator can take the decision of either 
going ahead with the integration process (using small or large steps); or going back- 
wards. Then, the coordinator computes the emission (using the numerical module, see 
section EPj) : and updates the current state of the automaton via the e variable (which 
is used to switch between two automaton states). 

The set of possible decisions (as shown in Tabled]), are based on two considerations: 

• In order to save computation time, we neglect the emission that does not con- 
tribute to the total brightness temperature. 

• On the other hand, we include, with a high spatial resolution, the emission of 
any structure in the solar atmosphere, which contribute to the total brightness 
temperature. 

2. The second component is the expert system which, based on the physical conditions of 
each specific point, decides if it is useful to integrate on this region and recommends 
the size of the following integration step. The recommendation is based on two plasma 
parameters: the plasma frequency [y v = 9x 10~ 3 y / n^ [MHz]) and the minimal emission. 

The plasma frequency is important to obtain the position (atmospheric height) of the 
interface between regions where electromagnetic waves, at any given frequency, can 
propagate or not. 

The minimal emission parameter defines the lower limit where the local emission is 
negligible and also controls both, the error due to this neglected flux; and the perfor- 
mance of the integration process, saving in this way a large amount of computation 
time. There is also another numerical error, associated to the small and large steps. 
We present the analysis of convergence of these errors in Appendix [Bj The minimal 
emission can be set by the user via "-min" parameter at the console or can be managed 
automatically by the code (see section . 

The expert system can recommends the integration steps (small or large), based on 
the following cases: 
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• If Up > is, then the wave can not propagate and the experts recommends a small 
integration step. We consider this case, because we want to know the height where 
the radiation starts propagating in the atmosphere. 

• If v p < v and the local emission is greater than the minimal emission (the amount 
of emission is important). The wave can propagate and the expert recommends 
small integration step (we want to analyze in detail the emission process). 

• \iv v <v and the local emission is lower than the minimal emission, the wave can 
propagate but there is not enough emission, therefore the expert recommends a 
large integration step (we want to save time in the computation process). 

The recommendations are managed by two variables: "q", the local behavior of the 
emission and "y", the size of the integration step (see Table [2]). These variables are 
transmitted to the coordinator, as well as the variable "state" which contains the 
current state of the Cellular Automaton (Tabled]). 

3. The Cellular Automaton is the logical structure that stores the stage of the integration 
process, is the memory of the system. It is controlled by two parameters: e (the variable 
that preserves the memory when the system switch from one state to another); and 
the stack, a logical structure which preserves a local memory of the automaton states. 
The stack help us to control the number of steps when the coordinator decides to "go 
back" in the integration process, and is represented by two integer variables: "i" and 
"n" (where n=largestep/smallstep, in this way, we warranty that the total length of 
the small and large steps is the same when the system enter in the "go back" process). 

The automaton can be in any of the four following states (see Table [3]): 

• Al: Integrating using small steps. 

• A2: Integrating using large steps. 

• A3: "Going back": I tried to integrate using large steps but I had to return 
because the local emission is larger than the minimal emission. Therefore I will 
integrate with small steps up to the returning point (this state shows the necessity 
of the Stack structure). 

• A4: Something is wrong. This is an error. 

In summary, the three components of Tulum conform a very efficient intelligent system 
to switch between integration steps; control the associated errors; and reproduce the emission 
with high spatial resolution. 
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3.2. 3D Geometrical Model 

The geometrical model was designed to optimize the computations of solar 3D struc- 
tures based on radial profiles of physical parameters (in general, quiet Sun models for the 
electronic density and temperature are given as radial profiles, starting at photospheric level 
and extending to different atmospheric altitudes). 

The origin of the coordinated system is located at the center of the solar sphere, the 
Z axis points towards the observer, the Y axis points to the solar North, and the X axis 
completes the system. In this geometry, a ray path from a given point in the plane of the 
sky to the observer is formed by a set of radial vectors (see Figure H]). These vectors describe 
both, the integration mesh and the radial values of density and temperature along this ray 
path. The radiative transfer equation is integrated along each ray path and the set of ray 
paths forms the 2D projection (on the plane XY) of the 3D model. 

In this geometry each point of the mesh is defined as: 

*o.a(z) = ( r (a x , Py, z), 6(ux, P y , z), (j)(a x , p y , z)) 

where r is the module of vector r ; 9 the angle between the Z axis and the projection of r*on 
the XZ plane; (f) is the angle between r and the ZY plane and z is the projection of r on to 
the Z axis. From the observer point of view, each ray path represents a pixel (x, y) on the 
projected 2D image and is defined by the angles a x and P y . Each ray path is divided in k 
points separated by a distance dl, for simplicity, we do not use directly dl but its projection, 
dz, on the Z axis. 

The mesh is defined by two constants: the Astronomical Unit, UA = 1.5 x 10 8 km and 
the solar radii R & = 6.96 x 10 5 km; plus the following variables: 

• n: The image resolution isnxn pixels. 

• x: Variable in the X direction ranging from — (n — l)/2 to (n — l)/2. 

• y: Variable in the Y direction ranging from — (n — l)/2 to (n — l)/2. 

• Rt- The maximum radial distance, in the 2D projection, considered for the integration 
(we use Rt = 2i? Q ). 

• F: Defines (in units of solar radii), the starting point of the integration process, F = 
means that the starting point lays in the plane XY; F = — 1 in a parallel plane located 
at one solar radii behind of the origin; and F = 1 in a parallel plane located at one 
solar radii in front of the origin (by default, we use F = —Rt)- 
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• H: is the final point of the integration process, the default value is H = Rt- 

• dl: integration step in km. 

The process starts by computing the matrix of angles: 

M n ,n = j(a(x),/%)) | <a:,y<^| 

Then, the initial and final integration points are computed for each element M XtV . These 
points are defined by the user (F and H, respectively). It is possible that some ray paths 
intersect the solar surface, for such cases we define F = zq, where zq is the projection of the 
intersection point on the Z axis. Once the initial and final integration points are known the 
code generates the set of points: 

L x , y = {r{a x , f3 y , z), 9(a x , (3 y , z), <fi(a x , f3 y , z) \ z < z < H and z = m*dz,m e N}, 

and solve Eq. [3j 



3.3. Model for thermal radio emission 



At quiet regions, the main contribution to the emission and absorption is due to free - 
free interactions, in particular, free electrons interacting with ions. The electron - electron; 
i on - ion; an d free - bound interactions, do not contribute significantly to the total emission 
(jDulki Il985). Even m ore, for radio emissions, only distant electron - ion i nteractions are 
important (]Dulklll985l ). Therefore, in this case, the absorption coefficient is (jDulki Il985l ): 



3c W v 2 m l / 2 {kTV 2 ) v/3 
where n« is computed using Saha equation ( Athay fc Thomas] 1961 ): 



(4) 



log H*±! = -0.1761 - log(P e ) + log + 2.5 logT - 

7k Ui 1 



(5) 



where Ui is the statistical weight; Xi is the ionization energy; P e = n e KT; and n e is the 
observed electronic density profile. 

Equation H] may be approximated, according to the appropriate Gaunt factor to: 

n P . fl8.2 + ln(T 3 / 2 )-ln^, T < 2 x ltf K 

> z, rii x < (6) 
1 24.5 + ln(T) - In v, T > 2 x 10 5 K 



Kjy w 9.78 x 10" 3 - 



v 



2^3/2 
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The source function is: 

2hu 3 1 

J = (7) 

c 2 exp{hu/kT) -1' y ] 

Although, at radio wavelengths hv << kT, it is possible to use the Rayleigh- Jeans approxi- 
mation: 

2ku 2 

h « — —T, 8 

c z 

Eq. [6] and |8] are solved by our model. We have simulated the solar emission, in the radio 
wavelengths range, and found a good agreement with observations (Section HJ). 



3.4. The Minimal Emission Parameter 

As shown in the upper panel of Figure El where we have plotted the total emission as 
a function of the photospheric height for different frequencies, from 7 GHz (black curve) 
to 7 THz (blue curve), above 3000 km the total emission have reached its final value for 
all frequencies. Obviously, this convergence occurs at different heights depending on the 
frequency, the minimum height of convergence (~ 590 km, marked with a vertical dotted 
line) corresponds to the 7 THz profile. We use this point as a reference height (h c ). 

On the central panel of Figure [3j we have plotted the "emission efficiency" , 

J eff = 1 - ex P(- r ^))> 

as a function of height for the same frequency range. Clearly the 7 THz profile has the 
lowest "emission efficiency" at all heights. Therefore, we can use the value of the "emission 
efficiency" of this profile at h c , this is, J e g = 1 x 10~ 4 , as the lower bound of the model 
(marked with a horizontal dotted line). The "emission efficiency" can reach lower values, at 
higher altitudes, but as seen in the upper panel, the contribution to the total emission (for 
the 7 THz profile) at these heights is negligible. 

As we do not know, before the computations, where the "emission efficiency" will reach 
this lower bound value, we check in the temperature model (thick line in the upper panel) 
and see that the temperature model reach its minimum value MIN(Tr) at h c . Therefore, by 
using Eq. [8] we are able to obtain the minimal significant emission, 

2ku 2 

J min < ^MIN(T^) x 10" 4 , (9) 
where MIN(Tr) is the minimum in the atmospheric temperature radial profile. 
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In order to show the correctness of the previous analysis, in the bottom panel of Figure 
[3] we have plotted the local emission profiles as a function of height for the same range of 
frequencies. The horizontal dotted line represents I min = 0A4K, computed using MIN(Tr) 
at 7 THz. As expected, this line intersects with 7 THz profile exactly at h c . 

As an example, we have marked (red line) a not negligible excess (i. e. above the dotted 
horizontal line) of local emission at 3 THz, from ~ 800 to ~ 1300 km of height (marked with 
vertical dashed lines). And, as shown in the upper panel by a red line in the 3 THz profile, 
only this excess contributes to the total emission. 

For lower frequencies we have marked with crosses the height where the local emission 
becomes negligible (this is, where each profile crosses the I min bound in the bottom panel). 
This height is also marked with crosses in the total emission profiles (upper panel), showing 
that the emission at each frequency already have converged to its final value, at the marked 
height . 

As shown by Figure H] where we have plotted the error associated with Eq. [9j for 
frequencies higher than ~ 40 GHz, the relative error of the final brightness temperature is 
lower than 1%. Whereas for lower frequencies the error is higher, due to the fact that there 
are large regions (at coronal heights) which contribute with low amounts of local emission. 



4. Results 



We compute the free-free thermal radio emission from an atmosphere of Hydrogen 
Helium gas, using published (radial) profiles of solar temperature and density. We per- 
formed a multi-frequency analysis, from 2 to 20 GHz, shown in Figure |5] by the continuous 



and lo ng-dash lines, and compared our results against observations reported by IZirin et al. 



( Il99ll ) (triangles) and similar published analysis. The continuous line is the output of our 
model using = n e in Eq. |H The long-dash line is the output of our mod el considering 



radiat ion from HII, Hell and Helll ions in Eq. |H The short-dash line is the iBastian et al. 



rli996h model which use sim i lar ph ysical considerations as our model. The dotted line is 



the lLandi fc Chiuderi Dragol ( 120031 ) model computed from the observ ed differenti al emission 
measure and using an empirical opacity function. We also plotted the lAllenl ( 119631 ) (dot-dash 
line) model. 

In Figure |6l we have plotted the brightness temperature difference, between observations 
and models, whit the same li ne code as Figure El This difference decreases with frequency. 
At ~ 5 GHz, our models and IBastian et al.l ( 119961 ) model have an excess of ~ 1.5 x 10 4 K, 
whereas lLandi &: Chiuderi Dragol ( 120031 ) model has an excess of ~ 2 x 10 4 K. At ~ 20 GHz, 
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all models have b etter agreeme n t with observations, although the excess computed by our 
models as well as iBastian et al.l ( 119961 ) model is only ~ 5 x 10 3 K. 



As our code uses a cellular automaton and an expert system to solve efficiently the 
radiative transfer equation, we are able to achieve integration times which are up to one 
order of magnitude shorter than direct integration codes (see Appendix [B]), this makes 
possible to generate high definition 2D images from 3D structures, in reasonably short times 
and using very short (1 km) integration steps. Therefore, Pakal can compute the emitting 
spectrum from highly detailed source structures, as the expected in new generation solar 
chromospheric models. 

We have c ompa red the performance of our code against a similar code published by 
Selhorst et al.l (120051 ) and against a linear integration process (See Appendix [B]) . We found 
that Pakal can improve the integration time up to one order of magnitude compared with 
the linear integration process and up to 1/3 when compared with lSelhorst et al.l (120051 ) code. 
We have performed a detailed analysis of the quiet Sun emission at 17 GHz simulated by 
P akal and using tempera ture and density profiles observed in UV and continuum (see details 



m 



de La Luz et al.l 120081 ). 



Figure [7] shows an equatorial cut of a 1024 by 1024 image of the computed quiet Sun 
emission at 43 GHz, where the limb brightening is clearly seen. In this case, we used 
integration steps of 10 km and a minimal local emission of 10~ 17 . We ran the code using 
the initial values shown in Appendix |A] The limb brightening show a maximum intensity of 
23000 K and a minimum of 8000 K. Observations at similar frequenc ies, made in the 1950's, 
repor ted brightness temperatures from 5700 K to 6000 K at 40 GHz dWhitehurst fc Mitchell 
19561 ). Although, based in later observations at 50 GHz (IReberl Il97ll ) predicted a higher 
emission at the center of the disk of 7500 K. 



5. Summary 

We have developed a new numerical code to solve the radiative transfer equation in a 
radial (3D) geometry for stellar atmospheres. The code is composed by four independent 
modules: i) numerical model; ii) geometry; iii) numerical methods; and iv) physical stellar 
models. This architecture allows easy changes when we want to test different physical models. 

We found that the minimum of the temperature profile, can be used to compute the 
lower boundary of the emission, this boundary guarantees the numerical convergence of the 
final brightness temperature. 
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By improving the geometry and the integration process, the code is able to reproduce, 
with better results, classical analysis of the solar radio emission, as the analysis of the 
depth of emission and multi-frequency analysis in ID; or 2D analysis of the limb brightening 
(Appendix Hj). 

The code is up to one order of magnitude faster than linear integration codes, and 
three times faster than similar published codes. As future work we are going to implement 
adaptative integration steps and develop the Message Passing Implementation (MPI) of the 
code which will work in multi processors-computers, with these improvements, the code will 
be able to solve the radiative transfer equation in non homogeneous structures with more 
complicated physical conditions as non-LTE, more chemical species and emission processes. 
Finally, the code is free under request to the author. 

This work was supported by UNAM-PAPIT grant INI 17309 and CONACyT grant 49395 
Thanks to Dr. R. Caballero for allow us to use his computer facilities. 
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A. Testing the model 

Pakal is able to deal with different opacity functions and chemical species with different 
ionization states. Although, in order to compare our code with previous published results, we 
ran the code twice, firstly, using a restricted initial value of densities, rij = n e (as shown by 
continuous lines in Figures [5] and [6]) , and secondly, considering the density of each ionization 
state of a diatomic gas formed by H-He (long-dash lines in Figures and E]). The code inputs 
where: 

1. Inputs from libraries: 

• Source function: Rayleigh- Jeans approximation (Eq. [S]). 

• Opacity: free- free emission (Eq. [6]). 

2. Inputs from files: 



This preprint was prepared with the A AS IATgX macros v5.2. 
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Radial profiles of temperature , elect ronic and Hydrogen densities: Here we use 
the model C of IVernazza et al.l (Il98ll ). for chrom ospheric and lo w transition zone 
heights . For c oronal heights, we use the model of iGabriell (119761 ) and reported by 



Foukall (11990( 1 . 



• Assuming He = 0.1 * H. 
3. Console inputs: These inputs changed for each particular simulation. 



B. Analysis of Convergence 

The convergence test is necessary when we want to prove the adequate functionality of 
any code. We have developed three convergence tests, which also help us to test the efficiency 
of the code. The first one involves the -detail parameter, which determines the length of 
the integration step when the code is performing a detailed integration process. The second 
test involves the -min parameter, which sets the minimum emission considered by the code 
(i. e., emission below this value is neglected). Finally, we analyze the -big parameter, which 
determines the length of the large integration step used when the local emission is negligible. 

In order to found a lower boundary for the minimum emission parameter (-min), we 
ran several simulations at different frequencies. We found that the -min parameter have not 
negligible effects when it is greater than the emission computed at 1% of the minimum of the 
temperature profile. Note that the final error of the model is associated to this parameter, 
at least the total error will be comparable to the minimum emission parameter and depends 
indirectly on the minimum step of integration. 

The combination of these parameters determine the efficiency of the code. If we use a 
very small number for the -detail parameter, the code will take long time for the integration 
process. On the other hand, the code will loose valuable information by using too large 
numbers in the -big or -min parameters. Therefore, we need to look for the best parameters 
in terms of the integration time and the stability of the output. In Figure Owe have plotted 
the brightness temperature (continuos lines) and the integration time (dotted lines), versus 
the varying parameter (-detail, -big and -min, respectively) so we can test the stability and 
performance of the code. The main idea is to find out the best parameters in of the shortest 
integrations times, but without affecting the final brightness temperature. Those tests were 
carried out by computing the emission over a ray path in a single pixel at position (0,0), i. 
e., in the center of the solar disk. 

The upper panel of Figure [8] shows the first test, the computed brightness temperature as 
a function of the small integration step (-detail parameter) using a constant large integration 
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step of 100 km. If we set the -detail parameter to 100 km, (i. e., the detail integration and 
the big integration steps are equal), the resultant algorithm is really poor, because it is 
integrating sequentially. In this case, the integration time (dotted line in the upper panel of 
Figure ED is very fast (< 1 sec), but the brightness temperature computed is far away from 
the right value (1.6 x 10 4 K). When the detailed integration steps are lower than 20 km, the 
emission converges to 1.6 x 10 4 i^, although the integration time grows exponentially. As 
instance, to generate an image of 1024 by 1024 using a small integration step of 10 km, the 
integration time is almost two months. If the small integration step is 1 km, the integration 
time will be around two years. 

To perform the second test, we left constant both: the small (0.5 km) and large (100 
km) integration steps, and allow variations of the minimal emission (-min). The continuos 
line in the middle panel of Figure El shows that the brightness temperature converges when 
the minimal emissions is lower than 10~ 13 . When the minimal emission is higher than 10~ 13 , 
the brightness temperature diverges and the integration times are shorter. As instance, if 
the minimal emission is 10~ 17 and the integration step of 0.5 km, an image of 1024 by 1024 
takes 85 days of integration. 

In order to find out the best value for both integration steps (third test), in the buttom 
panel of Figure El we have plotted the brightness temperature (continuos line) as a function 
of the large integration step , in terms of the small integration step (large = n x small), 
setting the minimal emission as 10 -17 and the small integration step as 0.5 km. 

Changes of the large integration step do not affect appreciably the brightness temper- 
ature, although, the integration time is largely affected by such changes. In this case, the 
integration time may vary in one order of magnitude. The minimum time of integration is 
reached at 60 km (see the dotted line of the button panel of Figure EJ), this is: 

&i<?[km] = 60 * detail = 30km. 

Performing the convergence analysis, but using the best parameters, is possible to obtain 
integration times which are one order of magnitude lower than direct integration process, as 
shown in Figure [9] For instance, the integration time for a 1024 by 1024 image with small 
integration steps of 10 km is now 11 days instead of two months. If the small integration 
step is 1 km Pakal now takes 39 days instead of two years (a super computer with 1024 
processors will take one hour to generate this image). 

We have c ompa red the performance of our code against a similar code published by 



Selhorst et al.1 (120051 ) which practically solves the same task but based on linear integration 



method. The border conditions are: 

• Spatial resolution: 1100 points = 770000 km. 
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• Steps of Integration: 800. 

• Size of a step: 50 km. 

The time results are as follows: 

• Selhorst code: 13 minutes. 

• Pakal code: 33 minutes. 

Changing the Pakal integration parameters to the optimal value, 

• Detail integration step: 50 km. 

• Large integration step: 100*50 km. 

• Minimal emission parameter: le-20. 

We obtain 

• Pakal code: 4 minutes. 



Summarizing we improved in one order of magnitude our results and in 1/3 the linear 
i nteg ration used in another method. Even more, that the method used by ISelhorst et al. 
( 120051 ) does not include the computation of the geometry of the problem. 
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Fig. 1. — General geometry used in the numerical model. 
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Automaton A2: Large Steps (Searching...) 




Fig. 2. — Detailed geometry used in Pakal showing, in an (Earth - P) ray path, three possible 
states of the integration process. 
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Fig. 3. — Multi-frequency spectrum from 7.5 GHz (black or upper curve) to 7 THz (blue 
or lower curve). The total emission (upper panel), "emission efficiency" (middle panel) and 
local emission (bottom panel) as a function of height from the photosphere. The thick black 
line in the upper panel is the radial temperature model. 
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Fig. 4. — Relative error in the final brightness temperature using min = I min parameter. 
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Fig. 5. — Quiet Sun brightness tem perature as a func tion of frequency from different models. 
Triangles show the observations of IZirin et al.l (1199 lh . 
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Fig. 6. — Brightness temperature difference between observations and models as a function 
of frequency. 
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Fig. 7. — A equatorial cut from 1024x1024 simulation image of the Sun at 43 GHz. The 
figure show the classical limb brightening. 
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Convergence Analysis 
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Fig. 8. — Analysis of convergence of Pakal. The brightness temperature (continuos lines) 
and the time of the integration process (dotted lines) for the variation of "-detail" (upper 
panel), "-min" (middle panel), and "-big" (bottom panel) parameters. 
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Fig. 9. — Analysis of convergence of Pakal. The brightness temperature as function of the 
short integration steps and using the optimal parameters; calling sequence: "./pakal -xy 
-nu 17e9 -min le-17 -detail X -big 60". The bottom panel shows the respective time 
(continuous line)and the time using linear integration (dotted line). 
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Table 1: Decision Table of the Coordinator. Where I is the local emission after the computa- 
tion; I is the incoming emission; S is the source function; r the local opacity; x a and Xb the 
two spatial coordinates; dzdetail the small integration step; and dzbig the large integration 
step. Using this decision table the coordinator chooses the integration step and computes 
the local emission /. 
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Variable 


Value 


Meaning 


y 





Next step small. 




1 


Next step large. 
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The wave can not propagate. 




1 


The emission is not enough. 




2 


There is enough emission. 



Table 2: State of the expert system, where q and y are the possible recommendations. 
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Table 3: Table of States of the Cellular Automaton. In this case n = dzbig / dzdetail. The 
variable "state" represents the memory stages of the integration process; "e" controls the 
switch between states; and the variable "i" is a stack into the automaton. 



